Dane pochodzą z Protein Data Bank, zawierają informację o ligandach. W analizie pominęliśmy kolumnę title, która zawierała powtórzone dane z innych kolumn oraz kolumnę weight_col, która składała się z wartości pustych. Dodatkowo w niektórych wierszach brakowało informacji odnośnie jednego z progów odcięcia. Predykcja ilości atomów oraz elektronów powiodła się z dość dużą trafnością co jest pokazanego w wynikach.
library(dplyr)
library(ggplot2)
library(ggforce)
library(gganimate) #devtools::install_github('thomasp85/gganimate')
library(tidyr)
library(caret)
library(DT)
library(summarytools)
library(gifski)
library(png)
library(data.table)
Aby zapewnić powtarzalność wyników ustawiamy stan losowego generatora liczb.
set.seed(123)
Do wczytania danych używamy funkci “fread”, aby zapewnić szybsze wczytanie danych. Kolumna title jest połączeniem kolumn pdb_code, res_name, res_id oraz chain_id, zatem możemy ją usunąć podczas wczytywania, aby uniknąć powtarzania informacji.
all_data <- fread("all_summary.csv", header = TRUE, dec=".", stringsAsFactors = FALSE) %>% select(-(blob_coverage:title))
Wiersze, które w kolumnie “res_name” zawierają niepożądaną przez nas wartość zostają usunięte.
cleaned_data <- all_data %>% filter(!res_name %in% c('UNK', 'UNX', 'UNL', 'DUM', 'N', 'BLOB', 'ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLN', 'GLU', 'GLY', 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'MSE', 'PHE', 'PRO', 'SEC', 'SER', 'THR', 'TRP', 'TYR', 'VAL', 'DA', 'DG', 'DT', 'DC', 'DU', 'A', 'G', 'T', 'C', 'U', 'HOH', 'H20', 'WAT'))
W zbiorze danych występuje kolumna, która równa jest NA we wszystkich wierszach.
rows_without_na_in_weight_co = filter(cleaned_data, !is.na(weight_col))
dim(rows_without_na_in_weight_co)[1]
## [1] 0
Jak widać brak wierszy, których kolumna weight_col nie ma wartości pustej, dlatego możemy ją wykluczyć z dalszej analizy. Usuwamy także wiersze zawierające wartość pustą.
cleaned_data_without_empty_col <- select(cleaned_data, -weight_col)
cleaned_data_without_empty_col <- cleaned_data_without_empty_col[complete.cases(cleaned_data_without_empty_col),]
Zebrane dane zawierają… 408 kolumn oraz 525666 wierszy. Dane są typu character, numeric, integer. Większość kolumn jest numeryczna. Ich podstawowe statystyki prezentują się tak:
Natomiast pozostałe kolumny są następujące:
| pdb_code | res_name | res_id | chain_id | skeleton_data | fo_col | fc_col | |
|---|---|---|---|---|---|---|---|
| Length:525666 | Length:525666 | Length:525666 | Length:525666 | Length:525666 | Length:525666 | Length:525666 | |
| Class :character | Class :character | Class :character | Class :character | Class :character | Class :character | Class :character | |
| Mode :character | Mode :character | Mode :character | Mode :character | Mode :character | Mode :character | Mode :character |
Naszą analizę ograniczymy do 50 najpopularniejszych wartości kolumny res_name.
get_top_n_res_names <- function(input_data, top_count) {
input_data %>%
select(res_name) %>%
group_by(res_name) %>%
count() %>%
arrange(desc(n)) %>%
head(top_count)
}
top_50_res_name <- get_top_n_res_names(cleaned_data_without_empty_col, 50)
data_with_most_common_res_names <- cleaned_data_without_empty_col %>% filter(res_name %in% top_50_res_name$res_name)
Rozkład jej wartości prezentuje się następująco
W celu sprawdzenia korelacji użyjemy korelacji Rho Spearmana, ponieważ rozkład wartości przynajmniej jednej kolumny nie jest rozkładem normalnym
correlation <- as.data.frame(as.table(cor(data_with_most_common_res_names %>% select_if(is.numeric), use="complete.obs", method="spearman")))
Usuniemy teraz korelacje kolumn samych ze sobą
correlation <- correlation %>%
rename(first_column = Var1, second_column = Var2, freq = Freq) %>%
filter(first_column != second_column)
Grupujemy po pierwszej kolumnie oraz dla każdej wartości obliczamy maksymalną wartośi. Następnie wyznaczmy 10 kolumn z największą korelacją oraz filtrujemy dane do wizualizacji
top_correlated <- correlation %>%
group_by(first_column) %>%
summarise(max=max(freq, na.rm = TRUE)) %>%
arrange(desc(max)) %>%
head(10)
correlation <- correlation %>% filter((first_column %in% top_correlated$first_column & second_column %in% top_correlated$first_column))
Wykres przedstawia ilości segmentów maski kształtu oraz maski gęstości elektronowej dla każdego progu odcięcia intensywności na podstawie danych z trzema najpopularniejszymi nazwami zasobu (res_name). Skala osi pionowej jest logarytmiczna, ponieważ różnice w wartościach są dość duże.
Kolumny do regresji zostaną wybrane na podstawie współczynnika korelacji Rho Spearmana. Następnie dane zostaną podzielone na zbiór treningowy oraz testowy w stosunku 70:30. Model zostanie oceniony na podstawie 5-krotnej, podwójnej oceny krzyżowej. Dodatkowo dokonamy predykcji na danych testowych oraz obliczymy miary R^2 oraz RMSE.
correlation <- cor(data_with_most_common_res_names %>% select_if(is.numeric), use="complete.obs", method="spearman")
correlation <- gather(as.data.frame(correlation)['local_res_atom_non_h_electron_sum',], 'column', 'correlation_ratio') %>%
filter(correlation_ratio > 0.8)
electron_count_predict_data <- select(data_with_most_common_res_names, correlation$column)
indexes <- createDataPartition(electron_count_predict_data$local_res_atom_non_h_electron_sum,p=0.7, list=F)
training_data <- electron_count_predict_data[indexes,]
testing_data <- electron_count_predict_data[-indexes,]
ctrl <- trainControl(
method = "repeatedcv",
number = 2,
repeats = 5)
fit <- train(local_res_atom_non_h_electron_sum ~ .,
data = training_data,
method = "lm",
trControl = ctrl)
fit
## Linear Regression
##
## 245042 samples
## 8 predictor
##
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times)
## Summary of sample sizes: 122520, 122522, 122521, 122521, 122521, 122521, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.46865 0.999972 0.1390262
##
## Tuning parameter 'intercept' was held constant at a value of TRUE
predicted_values <- predict(fit, newdata = testing_data)
rmse_electrons <- RMSE(testing_data$local_res_atom_non_h_electron_sum, predicted_values)
r2_electrons <- R2(testing_data$local_res_atom_non_h_electron_sum, predicted_values)
correlation <- cor(data_with_most_common_res_names %>% select_if(is.numeric), use="complete.obs", method="spearman")
correlation <- gather(as.data.frame(correlation)['local_res_atom_non_h_count',], 'column', 'correlation_ratio') %>%
filter(correlation_ratio > 0.8)
atom_count_predict_data <- select(data_with_most_common_res_names, correlation$column)
indexes <- createDataPartition(atom_count_predict_data$local_res_atom_non_h_count, p=0.7, list=F)
training_data <- atom_count_predict_data[indexes,]
testing_data <- atom_count_predict_data[-indexes,]
ctrl <- trainControl(
method = "repeatedcv",
number = 2,
repeats = 5)
fit <- train(local_res_atom_non_h_count ~ .,
data = training_data,
method = "lm",
trControl = ctrl)
fit
## Linear Regression
##
## 245042 samples
## 10 predictor
##
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times)
## Summary of sample sizes: 122521, 122521, 122522, 122520, 122521, 122521, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.03612161 0.9999924 0.00517815
##
## Tuning parameter 'intercept' was held constant at a value of TRUE
predicted_values <- predict(fit, newdata = testing_data)
rmse_atoms <- RMSE(testing_data$local_res_atom_non_h_count, predicted_values)
r2_atoms <- R2(testing_data$local_res_atom_non_h_count, predicted_values)
| Predykcja | RMSE | R^2 | |
|---|---|---|---|
| Ilości atomów | 0.0345164 | 0.9999931 | |
| Ilości elektronów | 0.5205734 | 0.9999658 |
top_3_res_name <- get_top_n_res_names(cleaned_data, 3)
res_name_classification_data <- data_with_most_common_res_names %>% filter(res_name %in% top_3_res_name$res_name)
res_name_classification_data <- select(res_name_classification_data, (part_00_shape_segments_count:part_02_density_Z_4_0), (resolution:FoFc_max), res_name) %>% mutate(res_name=factor(res_name))
indexes <- createDataPartition(res_name_classification_data$res_name,
p=0.7, list=F)
training_data <- res_name_classification_data[indexes,]
testing_data <- res_name_classification_data[-indexes,]
ctrl <- trainControl(
method = "repeatedcv",
number = 2,
repeats = 5)
fit <- train(res_name ~ .,
data = training_data,
method = "rf",
trControl = ctrl,
ntree=7)
fit
## Random Forest
##
## 84242 samples
## 324 predictor
## 3 classes: 'EDO', 'GOL', 'SO4'
##
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times)
## Summary of sample sizes: 42122, 42120, 42121, 42121, 42120, 42122, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 2 0.6729873 0.4886008
## 163 0.6799530 0.5000081
## 324 0.6795162 0.4993857
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 163.
predicted_classes <- predict(fit, newdata = testing_data)
confusionMatrix(data = predicted_classes, testing_data$res_name)
## Confusion Matrix and Statistics
##
## Reference
## Prediction EDO GOL SO4
## EDO 4445 2868 1030
## GOL 2909 6544 1667
## SO4 1043 1772 13824
##
## Overall Statistics
##
## Accuracy : 0.6873
## 95% CI : (0.6825, 0.6921)
## No Information Rate : 0.4576
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.5113
## Mcnemar's Test P-Value : 0.3107
##
## Statistics by Class:
##
## Class: EDO Class: GOL Class: SO4
## Sensitivity 0.5294 0.5851 0.8368
## Specificity 0.8593 0.8164 0.8562
## Pos Pred Value 0.5328 0.5885 0.8308
## Neg Pred Value 0.8576 0.8143 0.8614
## Prevalence 0.2326 0.3098 0.4576
## Detection Rate 0.1231 0.1813 0.3829
## Detection Prevalence 0.2311 0.3080 0.4609
## Balanced Accuracy 0.6943 0.7007 0.8465